Thermal properties of fluorinated graphene 
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Large scale atomistic simulations using the reactive force field approach (ReaxFF) are imple- 
mented to investigate the thermomechanical properties of fluorinated graphene (FG). A new set 
of parameters for the reactive force field potential (ReaxFF) optimized to reproduce key quantum 
mechanical properties of relevant carbon-fluor cluster systems are presented. Molecular dynamics 
(MD) simulations are used to investigate the thermal rippling behavior of FG and its mechanical 
properties and compare them with graphene (GE), graphane (GA) and a sheet of BN. The mean 
square value of the height fluctuations (h 2 ) and the height-height correlation function H{q) for 
different system sizes and temperatures show that FG is an un-rippled system in contrast to the 
thermal rippling behavior of graphene (GE). The effective Young's modulus of a flake of fluorinated 
graphene is obtained to be 273 N/m and 250 N/m for a flake of FG under uniaxial strain along 
arm-chair and zig-zag direction, respectively. 



I. INTRODUCTION 

The fascinating properties of single layer graphene 
(GE) have triggered a broad interest in the solid state 
physics community [1—5]. Despite its high electron mo- 
bility [6] , the zero band gap defies its employment in nano 
transistors where it is desirable to have a large on-off ra- 
tio between conducting and non-conducting states. A 
band gap can be induced by the addition of adatoms, 
which changes locally the hybridization of the carbon 
(C) atoms, but also modifies the electron mean free-path 
affecting the electron transport properties. Hydrogen 
(H) and fluor (F) are two well tested candidates [7-10] 
that leads to a large band gap opening. Graphane (GA, 
hydrogenated graphene) and fluorographene (FG) have 
been studied both experimentally and theoretically [11- 
14] to engineer the band gap. 

When H or F atoms are attached to the C atoms of 
GE, the bonds transit from an sp 2 to an sp 3 hybridiza- 
tion, which turns the conjugated graphitic C-C bonds 
into single C-C bonds. In the fully covered cases both GA 
and FG are insulating materials [7, 8] and the structure 
changes locally the planar shape of GE into an angstrom 
scale out-of-plane buckled shaped membrane [15] known 
as chair configuration [16, 17]. 

From its potential applications in nanotechnology, FG 
is a more favourable material than GA since the C- 
F bonds are energetically more stable than the C-H 
bonds [13, 15, 17, 18]. Fluorographene has a very large 
temperature-dependent resistance and when the fluor 
content is increased it induces large changes in the elec- 
tron transport [19]. As in GE, it is expected that temper- 
ature also affects strongly the lattice structure and the 
mechanical properties of FG. 

According to the Mermin- Wagner theorem, [20], ther- 



mal excited ripples in two dimensional-like materials 
(GE, bilayer GE, GA and FG) has to play an impor- 
tant role in the stability of the membrane. While in 
GE and bilayer GE the corrugations are well described 
within the theory of two dimensional continuous mem- 
branes [21, 22], for GA instead, we recently found that 
the angstrom scale buckling of the carbon layer of GA 
prevents the formation of intrinsic long wavelength ther- 
mal ripples for temperatures up to at least 900 K [23]. 

Since the C atom has a higher (lower) electronegativity 
than H (F), it will take (give) away charge from the H 
(F) atom and consequently transforms the resulting C-H 
and C-F covalent bonds into sp 3 bonds. Therefore, it is 
expected that similar rippling effects as in GA will occur 
in FG although the C-F bonds are somewhat stronger 
than the C-H bonds. The latter is due to the larger 
amount of charge that is shifted from C to F as com- 
pared to the one from H to C [13]. However in order 
to simulate large scale FG samples an appropriate force 
field is needed which describes the true chemical bond 
in C-F. Indeed, the absent of such a suitable interatomic 
potential for C-F restricted most of the recent studies to 
ab-initio calculations of their electronic properties using 
a small computational unit cell. 

ReaxFF potential serves to describe both bond and 
non-bond interactions in solids. Recently, such poten- 
tials were parameterized and were well tested for dif- 
ferent kind of structures, e.g. hydrocarbons [24], car- 
bon allotropes [25], etc. In this study we present a new 
set of parameters for ReaxFF, appropriate for structures 
with C-F bonds. Using molecular dynamics (MD) sim- 
ulations over large scale samples we study the thermal 
corrugations of FG and compare the results with those 
found for GA, GE and BN. We show that fully covered 
FG follows the same trend as GA and does not develop 



long-wavelength ripples or significant corrugation. The 
bending rigidity n of FG is found to be larger than the 
one of GE, GA and BN. Furthermore, k turns out to be 
temperature independent. Our results indicate that long- 
wavelength ripples are instead present in partial covered 
FG samples with a larger amplitude as compared to GA. 

The paper is organized as follows. In Sec. II, we intro- 
duce a new set of parameters for the ReaxFF potential 
of the C-F covalent bond. Then, in Sec. Ill using the 
introduced parameters, we analyze the thermal rippling 
behavior. Here, we consider both fully and partially cov- 
ered graphene sheets by F atoms. All the results are 
compared with those previously found for graphane. We 
also estimate the effective Young's modulus of FG flakes. 
We conclude the paper in Sec. IV. 
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II. REAXFF POTENTIAL FOR 
FLUOROGRAPHENE 

ReaxFF [24] is a general bond-order dependent poten- 
tial that uses a relationship between bond distance and 
bond order on the one hand and a relationship between 
bond order and bond energy on the other hand to de- 
scribe bond formation and dissociation. Many body in- 
teractions such as the valence angle and torsional inter- 
actions are formulated as function of bond order so that 
their energy contributions vanish smoothly upon bond 
dissociation. Non bonded interactions, namely Coulomb 
and van der Waals interactions, are calculated between 
every pair of atoms irrespective of their connectivity. Ex- 
cessively close range interactions are avoided by shield- 
ing. ReaxFF uses the geometry dependent charge calcu- 
lation scheme (EEM scheme) of Mortier et al [26] . The 
system energy in ReaxFF consists of a sum of terms: 

Esys — ^bond H~ ^under H~ ^over H~ H~ ^val H~ -^penH - 
-Etors H~~ E con j -\- E vc i\y aa i s -\- E(j ou i omj \). 

A detailed description of each of these terms and 
their functional forms can be found in the original 
work [24]. The reactive force field for C/F contain- 
ing systems was developed by parameterizing the poten- 
tial against DFT data obtained at the B3LYP/6-31g** 
level of theory (which is implemented in Schrodinger [27] 
which is an electronic structure packages) for various 
quantities such as fluorine and carbon atom charges 
in H 3 C-CF 2 -CH 3 , C-F and C-C bond lengths in 
H 3 C-CF 2 -CH 3 and H 3 C-CF(CH 3 )-CH 3 , F-F bond 
length in the F 2 molecule, potential energy curve for 
C-F bond dissociation in H 3 C-CF 2 -CH 3 , F-C-F an- 
gle bending in H 3 C— CF 2 — CH 3 , C-C-F angle bending 
in H 3 C-CF 2 -CH 3 and F-C-C-F dihedral twisting in 
F 2 C=CF 2 along with various chemical reactions involv- 
ing fluoroalkanes and fluoroalkenes. The results of the 
force field training are presented in Figs. l(a)-(d) and in 
Table II. Fig. 2 depicts the geometrical quantities rele- 
vant to Figs. l(a)-(d). It can be seen from Table II that 
ReaxFF reproduces closely the DFT based equilibrium 



FIG. 1: (Color online) Comparison between DFT (solid 
squares) and ReaxFF (open circles) results for: (a) C-F bond 
dissociation in H3C-CF2-CH3 (b) C-C-F angle bending in 
H 3 C-CF(CH 3 )-CH 3 , (c) F-C-F angle bending in H 3 C-CF 2 - 
CH 3 ,and (d) F-C-C-F dihedral twisting in F 2 C=CF 2 . 



geometries for various compounds. ReaxFF predicts F 2 
dissociation energy of 36.6 kcal/mol, in very good agree- 
ment with the DFT value of 37 kcal/mol. Fig. 1(a) shows 
that ReaxFF based potential energy curve for the C-F 
bond dissociation in H 3 C— CF 2 — CH 3 closely follows the 
DFT based potential energy curve. ReaxFF is able to 
predict very precisely the equilibrium C-F bond length 
(see Table II) and the C-F bond dissociation energy. Sim- 
ilarly the force field can closely reproduce the DFT based 
potential energy curve and the equilibrium geometry (see 
Table II) for C-C-F angle bending and the C-F-C angle 
bending as shown in Figs. 1(b)- (c). Figure 1(d) shows 
the variation of the potential energy upon F-C-C-F dihe- 
dral angle twisting. Though ReaxFF predicts the correct 
trend, the torsional rotation barrier in ReaxFF is around 
18 kcal/mol lower than that predicted by DFT. Over- 
all, the ReaxFF force field for C/F systems can closely 
reproduce the DFT based energies and geometries for a 
number of molecules and reactions. This force field will 
now be employed in large scale fully reactive molecular 
dynamics simulation of C/F containing systems. 



In the next section we study the thermal structural 
fluctuations and mechanical properties of a single layer of 
FG using large scale atomistic simulations employing the 
presented ReaxFF parameters. These parameters were 
implemented in the large-scale atomic/molecular mas- 
sively parallel simulator package LAMMPS [28, 29]. 
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TABLE I: Comparison of equilibrium geometrical parameters between ReaxFF and DFT. 

DFT ReaxFF 

F 2 bond length 1.43 A 1.4012 A 

C-F bond length in H 3 C-CF 2 -CH 3 1.384lA 1.4057A 
C-F bond length in H 3 C-CF(CH 3 )-CH 3 1.384lA 1.4158A 
Non-bonding C-F distance in CF 2 dimer 2.00 A 2.4471A 
F-C-F angle in H 3 C-CF 2 -CH 3 105.65° 107.2197 ° 

C-C-F angle in H 3 C-CF(CH 3 )-CH 3 106.2 ° 109.9625 ° 



(a) 



(b) 




(c) 




(d) 




FIG. 2: (Color online) (Color online) The molecules used for 
parameterizing the ReaxFF force field in this study, (a) C-F 
bond in H 3 C-CF 2 -CH 3 (b) F-C-F angle in H 3 C-CF 2 -CH 3 (c) 
C-C-F angle in H 3 C-CF(CH 3 )-CH 3 , and (d) F-C-C-F dihedral 
angle in F 2 C=CF 2 . The atoms constituting the geometric 
parameters are represented by balls while the rest of the atoms 
are represented by sticks. F atoms are colored brown, C atoms 
are colored green and H atoms are colored blue. 



III. RESULTS 
A. Thermal rippling behavior of FG 

In order to study the rippling behavior of FG we con- 
sidered a square shaped computational unit cell of FG 
with both armchair and zigzag edges in the x and y di- 
rections. Partial fluor contents of 10%, 50%, 70%, 90% 
and the fully covered 100% case (with a total number of 
N = 17280 atoms) were studied. In our simulation we 
employed the NPT ensemble with P=0 using the Nose- 
Hoover thermostat and varied the temperature from 10 K 
to 900 K. Figure 3 shows the obtained buckled shape of 
fully fluorinated sample after relaxation which is in agree- 
ment with recent DFT results [13]. 

One would expect that the thermal excited ripples in 
FG can be described by membrane theory for a 2D con- 
tinuous membrane [30] . This theory, described in a series 
of related works [23, 31-33], is supposed to be universal 
and independent of the atomic scale details of the mem- 
brane. The main predictions of this theory are as fol- 




FIG. 3: (Color online) (a) Side view of the buckled structure, 
known as chair configuration, of fully fluorinated graphene. 
The averaged bond-angle, and C-C and C-F distances, are 
respectively 109.5°, dc-c =1.58 A and dc-F =1.41 A at 
room temperature. 



lows. Let h be the out-of plane displacement of a given 
atom of a sheet, then the Fourier transform of the height- 
height correlation function is in the harmonic approxima- 
tion given by 



H(q) = <|%)| 2 



Nk B T 
K,S q 4 



(1) 



where q is the wave- vector, N is the number of atoms, 
So is the surface area per atom, n is the bending rigidity 
of the membrane, and ks is the Boltzmann constant. 

In the large wavelength limit, anharmonic couplings 
between bending and stretching modes are important re- 
sulting in a renormalization of the (^-dependent behavior 



H(q) 



Nk B T 



(2) 



where 77 is an universal scaling exponent which is about 
« 0.85 [34-36]. 

In order to compare our results for FG with other 
two dimensional materials, we used a modified Tersoff 
potential (which is an ordinary defined potential in the 
LAMMPS package [29]) according to the parameters pro- 
posed by Sevik et al for the h-BN sheet [37]. To simulate 
GE and GA we have used the AIREBO potential [38] 
which is suitable for simulating hydrocarbons. 

Recently, we found that in GA, H(q) acquires a strong 
renormalization for small wave-vectors q and the layer 
remains almost flat even for temperatures as high as 900 
K [23]. Here we will analyze the thermal rippling be- 
havior of FG and compare it with GA. A comparison 
with GE and BN single layers which behave as 2D mem- 
branes [31, 39] will also be presented. H(q) for FG was 





AIREBO 


ReaxFF 




n(eV) {h 2 )(A 2 ) 


K (eV) (h 2 )(A 2 ) 


GE 


1.165 1.307 


1.16 0.627 


GA 


10.19 0.070 


7.26 0.074 


GF 




5.07 0.073 



TABLE II: 

Comparison of AIREBO and ReaxFF for the bending 
rigidity and (h 2 ) for GE, GA and FG at 300 K. 

calculated following the steps described in our previous 
work [23]. 

Starting from a pure GE sheet, the variation of the 
height-height correlation function H(q) at room temper- 
ature for different partial fluor contents is shown in Fig. 4 
(a). The curves were shifted for a better comparison. We 
found that while for 10 to 90% coverage, H(q) follows 
Eq. (2) up to small q— values which is similar to the case 
of GE [32]. But, for fully FG at q* w 0.2 A~\ H(q) de- 
viates from the harmonic law (solid line) and approaches 
roughly a constant value similar to was previously found 
for GA [23]. In the inset of Fig. 4(a) we show the square 
average of the out-of-plane fluctuations (h 2 ) at 300 K. 
Notice that the out of plane fluctuations for partially 
covered samples are considerably larger for FG than for 
GA. The temperature dependence of H(q) for fully flu- 
orinated graphene is shown in Fig. 4(b). Irrespective of 
temperature, the short wave-length limit of H(q) tends 
always approximately to a constant value. The charac- 
teristic g- value where H(q) deviates from the harmonic 
approximation result decreases with increasing tempera- 
ture. 

The renormalization of H(q) for long wavelengths in- 
dicates the suppression of large out-of-plane height fluc- 
tuations. In Fig. 5(a) we compare the behavior of (h 2 ) 
against temperature for GE, BN (panel (a.l)), FG and 
GA (panel (a. 2)). Notice that (h 2 ) increases from 0.7 A 2 
up to 4 A 2 in BN and from 0.7 A 2 up to 2 A 2 in GE 
when temperature is varied from 10 K up to 900 K. Due 
to the absence of long wave-length ripples, (h 2 ) remains 
approximately constant for GA and FG, and the varia- 
tions are smaller than those for BN and GE, over the 
same temperature range. The temperature dependence 
of the bending rigidity ft, computed from the harmonic 
part of H(q) is shown in Fig. 4(b). Note that the larger 
magnitude for GA and FG is a consequence of the smaller 
corrugations present in these materials. We also find the 
opposite temperature dependence for BN and GE when 
compared with GA and FG. In this sense GE and BN 
behave anomalously. The corresponding bending rigid- 
ity and (h 2 ) at room temperature for GE, GA and FG 
are listed in Table II. 

Density functional calculations for band structure of 
FG (and GA) [40] shows that the acoustic out-of-plane 
modes (ZA) in FG (and GA) are different from that of 
GE. The most important difference from GE is the de- 
coupled optical and acoustic bands in FG and GA close 
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FIG. 4: (Color online) (a) Log-log plot of H(q) for different 
coverage of F atoms at T=300 K. The solid lines show the 
harmonic q~ 4 behavior valid in the limit of large q values. 
Note the strong deviation, starting roughly at q* w 0.2A -1 in 
the limit of long wave-lengths, for the case of fully fluorinated 
graphene. The variation of (h 2 ) is shown in the inset, (b) 
H(q) for fully fluorinated FG at different temperatures. 



to the T point. The light H atom contributes to the high- 
est phonon frequencies which is not the case for F atoms. 
It is also seen that the ZA modes close to the T point 
for FG (and GA) are not well fitted by a quadratic func- 
tion in contrast to the GE case. This is clear indication 
of anharmonicity. In summary, the atomistic details of 
the structure of FG is more complicated and therefore 
more details of this structure should be included in any 
continuum theory. 

Before ending this section, note that as we discussed 
in our previous work [23], the scaling with the system 
size present in GE is no longer valid for FG and GA 
((h 2 ) in FG and GA is almost constant irrespective to 
the system size). The lower wavelengths (q) adopted for 
the calculation of H(q) are equal to q x -min = j 1 and 
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FIG. 5: (Color online) Variation of (a) (/i 2 ), and (b) k against 
temperature for FG (open squares), GA (filled squared), GE 
(open circles), and BN (stars). 



Qy-min = and represent the 'computational' cut-offs 
of possible large wavelength ripples where l x and l y are 
the dimension of the system. Notice that deviation from 
the harmonic behavior takes place at larger value of q 
and hence this effect can not be a finite size effect and is 
instead an intrinsic phenomenon of the material. 



B. Effective Young's modulus 



In order to study the mechanical stiffness we consider 
an FG flake with dimension l x x l y = 15 x 15 nm 2 . Be- 
fore the stretching process, the sample is equilibrated for 
5 ps (i.e. 50,000 time steps). Stretching direction is al- 
ways along x and the uniaxial strain is applied within the 
NPT ensemble [41] where the pressure is slowly increased, 
i.e. 2GPas/ps. In this section the lateral edges (in the 
y direction) were taken as the arm-chair direction hav- 
ing both free (FBC) and periodic boundary conditions 
(PBC). We kept temperature fixed at T=10K. 

The total strain energy per atom of the strained flake 
can be written as a function of the imposed strain (e) 



So 
lit 



7(e) 



So 



Ye 2 



(3) 



where Eq is the energy of the infinite planar undeformed 
flake, 7(e) is the excess edge energy, and Y is Young's 
modulus (Y) of the flake. 

For nano-ribbons with no lateral edges we have 7 = 
(assuming that the longitudinal edges which are fixed 
make no contribution). This is due to the fact that free 
edges increase the energy due to buckling and bond-order 
changes [42]. Recently Lu et al used the Brenner poten- 
tial [43] in molecular dynamics simulations and studied 
the excess edge energy of graphene nano-ribbons as a 



function of width and chirality [42] . Our systems are dif- 
ferent from those of Ref. [42]. In contrast to Ref. [42] 
we are not interested in effects due to the edge energy 
effect and the size dependence. We rewrite Eq. (3) in the 
following as an effective Young's modulus which qualita- 
tively gives a good description of the mechanical stiffness 
of all the examined 2D materials. Nevertheless our re- 
sults are in qualitative agreement with those reported by 
Lu at al, i.e. increasing of total energy for the FBC case 
as compared to a nano-ribbon. Assuming a quadratic 
relation for 7(e) = l -^e 2 valid for small e, the simplest 
method to estimate Young's modulus is by fitting the 
quadratic function to the total energy (per area): 



Ej 



1 



Ye 



efft 



(4) 



where Y e ff is the effective Young's modulus of the sys- 
tem. Using aforementioned fitting process we found Y e f f 
respectively for a flake with arm-chair and zig-zag FG, to 
be 273 N/m and 250 N/m. Notice that the experimental 
result is 100 N/m for not perfect FG [8] while the DFT 
result is 250 N/m [18]. The latter disagreement between 
theory and experiment may be explained due to the fact 
that in experimental samples the fluor-to-carbon ratio is 
larger than unity, i.e. 1.1 [8], because of the presence of 
defects. Such defects become active regions which can 
adsorb the free F (and even H) atoms. Therefore, in 
the defected parts more F atoms will be found which is 
responsible for a F/C ratio larger than one. 

In order to understand the effect of the different bound- 
ary conditions, we depict in Fig. 6 the variation of Et 
per atom with e for flakes with both FBC (dashed lines) 
and PBC, i.e. nano-ribbon (solid lines). It is seen that 
for flakes with FBC the free edges result in an increase 
of the energy. The inset shows the difference between 
two curves, i.e. AEt = Efbc — Epbc which is posi- 
tive. Because the free boundaries have many dangling 
bonds which are not saturated by F atoms it results in 
extra energy. This can also occur in other systems, e.g. 
graphene [44]. Notice that for the studied low tempera- 
ture here, i.e. T=10 K we do not expect that bond recon- 
struction at the edges is important. Notice that even by 
saturating all the bonds by F, still the change in the bond 
order term in ReaxFF (due to different chemical environ- 
ment of the boundary atoms) results in higher energy as 
compared to PBC. 

Furthermore, both FBC and PBC results exhibit a 
quadratic behavior which is an indication of the harmonic 
regime for the applied uniaxial stain. As is clear from the 
inset of Fig. 6, the difference between the two curves is 
increasing with applied strain. This is due to the devia- 
tion from equilibrium for the C-F bonds, C-C-F (F-C-F) 
bond angles, and the dihedral angles (F-C-C-F torsion 
angle) of the free edge atoms. The larger the strain (and 
the larger the length of ribbon), the larger the deviation 
from equilibrium for the bonds and the angles. In the 
PBC case there is no such edge effect but nevertheless 
because of the fixing of the two longitudinal ends (the 
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FIG. 6: (Color online) (a) Variation of total energy against 
uniaxial strain for FG subjected to free boundary condition 
(FBC) and periodic boundary condition (PBC) for the lateral 
edges, i.e. the dashed and solid curves respectively. The inset 
shows the difference between the two energy curves. 



edges which are under uniaxial stress) the energy varia- 
tion of the PBC system should be different from that of 
an infinite FG which is periodic in both directions while it 



is under tension from the arm-chair direction. The fixed 
longitudinal ends do not have any effect in our results be- 
cause both FBC and PBC have the same contributions. 



IV. CONCLUSIONS 

We provided a new set of parameters for the ReaxFF 
potential for the C-F covalent bond and tested it on var- 
ious molecules. Subsequently, molecular dynamics sim- 
ulations were used to investigate the thermal rippling 
behavior and the mechanical response of fluorographene 
(FG) under uniaxial stress. The obtained results are 
compared with those for graphene (GE), graphane (GA) 
and hexagonal boron nitride sheet (BN). We found that 
fluorographene remains a flat sheet similar to graphane 
even at high temperature, i.e. up to 900 K. The bending 
rigidity of FG is found to be independent of tempera- 
ture and its Young's modulus is in good agreement with 
experiment. 
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